Basic fit setup and numerics

7. Basic fit setup and numerics#

7.1. Init and pulling data#

Here the setup is mainly handled by some basic scripts, these follow the outline in the PEMtk documentation [11], see in particular the intro to fitting.

import wget
from pathlib import Path

dataName = 'n2fitting'

# Pull N2 data from ePSproc Github repo
# URLs for test ePSproc datasets - n2
# For more datasets use ePSdata, see https://epsproc.readthedocs.io/en/dev/demos/ePSdata_download_demo_300720.html
urls = {'n2PU':"https://github.com/phockett/ePSproc/blob/master/data/photoionization/n2_multiorb/n2_1pu_0.1-50.1eV_A2.inp.out",
        'n2SU':"https://github.com/phockett/ePSproc/blob/master/data/photoionization/n2_multiorb/n2_3sg_0.1-50.1eV_A2.inp.out",
        'n2ADMs':"https://github.com/phockett/ePSproc/blob/master/data/alignment/N2_ADM_VM_290816.mat",
        'demoScript':"https://github.com/phockett/PEMtk/blob/master/demos/fitting/setup_fit_demo.py"}



# Set data dir
dataPath = Path(Path.cwd(), dataName)

# Create and pull files if dir not present (NOTE doesn't check per file here)
if not dataPath.is_dir():
    dataPath.mkdir()

    # Pull files with wget
    for k,v in urls.items():
        wget.download(v+'?raw=true',out=dataPath.as_posix())  # For Github add '?raw=true' to URL


# List files
list(dataPath.glob('*.out')) + list(dataPath.glob('*.mat'))

[PosixPath('/home/jovyan/jake-home/buildTmp/_latest_build/html/doc-source/part2/n2fitting/n2_1pu_0.1-50.1eV_A2.inp.out'),
 PosixPath('/home/jovyan/jake-home/buildTmp/_latest_build/html/doc-source/part2/n2fitting/n2_3sg_0.1-50.1eV_A2.inp.out'),
 PosixPath('/home/jovyan/jake-home/buildTmp/_latest_build/html/doc-source/part2/n2fitting/N2_ADM_VM_290816.mat')]

7.3. Compute AF-\(\beta_{LM}\) and simulate data#

With all the components set, some observables can be calculated. For testing, we’ll also use this to simulate an experiemental trace…

Here we’ll use self.afblmMatEfit(), which is also the main fitting routine, and essentially wraps epsproc.afblmXprod() to compute AF-\(\beta_{LM}\)s (for more details, see the ePSproc method development docs).

If called without reference data, the method returns computed AF-\(\beta_{LM}\)s based on the input subsets already created, and also a set of (product) basis functions generated - these can be examined to get a feel for the sensitivity of the geometric part of the problem, and will also be used in fitting to limit repetitive computation.

7.3.1. Compute AF-\(\beta_{LM}\)s#

# data.afblmMatEfit(data = None)  # OK
BetaNormX, basis = data.afblmMatEfit()  # OK, uses default polarizations & ADMs as set in data['subset']
# BetaNormX, basis = data.afblmMatEfit(ADM = data.data['subset']['ADM'])  # OK, but currently using default polarizations
# BetaNormX, basis = data.afblmMatEfit(ADM = data.data['subset']['ADM'], pol = data.data['pol']['pol'].sel(Labels=['x']))
# BetaNormX, basis = data.afblmMatEfit(ADM = data.data['subset']['ADM'], pol = data.data['pol']['pol'].sel(Labels=['x','y']))  # This fails for a single label...?
# BetaNormX, basis = data.afblmMatEfit(RX=data.data['pol']['pol'])  # This currently fails, need to check for consistency in ep.sphCalc.WDcalc()
                                                                    # - looks like set values and inputs are not consistent in this case? Not passing angs correctly, or overriding?
                                                                    # - See also recently-added sfError flag, which may cause additional problems.

7.3.2. AF-\(\beta_{LM}\)s#

The returned objects contain the \(\beta_{LM}\) parameters as an Xarray…

# Line-plot with Xarray/Matplotlib
# Note there is no filtering here, so this includes some invalid and null terms
BetaNormX.sel(Labels='A').real.squeeze().plot.line(x='t');

../_images/a36c9e6d410e51092c05ca332eec08a59cc2739357abf95148eaf2c2262206d0.png

… and the basis sets as a dictionary.

basis.keys()

dict_keys(['BLMtableResort', 'polProd', 'phaseConvention', 'BLMRenorm'])

7.4. Fitting the data#

In order to fit data, and extract matrix elements from an experimental case, we’ll use the lmfit library. This wraps core Scipy fitting routines with additional objects and methods, and is further wrapped for this specific class of problems in pemtkFit class we’re using here.

7.4.1. Set the data to fit#

Here we’ll use the values calculated above as our test data. This currently needs to be set as self.data['subset']['AFBLM'] for fitting.

# data.data['subset']['AFBLM'] = BetaNormX  # Set manually

data.setData('sim', BetaNormX)  # Set simulated data to master structure as "sim"
data.setSubset('sim','AFBLM')   # Set to 'subset' to use for fitting.

Subselected from dataset 'sim', dataType 'AFBLM': 156 from 156 points (100.00%)

# Set basis functions
data.basis = basis

7.4.2. Adding noise#

# Add noise with np.random.normal
# https://numpy.org/doc/stable/reference/random/generated/numpy.random.normal.html
# data.data['subset']['AFBLM']

import numpy as np
mu, sigma = 0, 0.05  # Up to approx 10% noise (+/- 0.05)
# creating a noise with the same dimension as the dataset (2,2) 
noise = np.random.normal(mu, sigma, [data.data['subset']['AFBLM'].t.size, data.data['subset']['AFBLM'].l.size])
# data.BLMfitPlot()

# Set noise in Xarray & scale by l
import xarray as xr
noiseXR = xr.ones_like(data.data['subset']['AFBLM']) * noise
# data.data['subset']['AFBLM']['noise'] = ((data.data['subset']['AFBLM'].t, data.data['subset']['AFBLM'].l), noise)
# xr.where(noiseXR.l>0, noiseXR/noiseXR.l, noiseXR)
noiseXR = noiseXR.where(noiseXR.l<1, noiseXR/(noiseXR.l))  # Scale by L

data.data['subset']['AFBLM'] = data.data['subset']['AFBLM'] + noiseXR
data.data['subset']['AFBLM'] = data.data['subset']['AFBLM'].where(data.data['subset']['AFBLM'].m == 0, 0)

data.BLMfitPlot()

Dataset: subset, AFBLM

../_images/8ce29b865817b39b9414464b5c56bfc93648c4ca8da96312bffa49a9f31b4229.png

7.4.3. Setting up the fit parameters#

In this case, we can work from the existing matrix elements to speed up parameter creation, although in practice this may need to be approached ab initio - nonetheless, the method will be the same, and the ab initio case detailed later.

# Input set, as defined earlier
# data.data['subset']['matE'].pd

# Set matrix elements from ab initio results
data.setMatEFit(data.data['subset']['matE'])

Set 6 complex matrix elements to 12 fitting params, see self.params for details.
Auto-setting parameters.

name	value	initial value	min	max	vary	expression
m_PU_SG_PU_1_n1_1_1	1.78461575	1.7846157536101068	1.0000e-04	5.00000000	True
m_PU_SG_PU_1_1_n1_1	1.78461575	1.7846157536101068	1.0000e-04	5.00000000	False	m_PU_SG_PU_1_n1_1_1
m_PU_SG_PU_3_n1_1_1	0.80290495	0.802904951323892	1.0000e-04	5.00000000	True
m_PU_SG_PU_3_1_n1_1	0.80290495	0.802904951323892	1.0000e-04	5.00000000	False	m_PU_SG_PU_3_n1_1_1
m_SU_SG_SU_1_0_0_1	2.68606212	2.686062120382649	1.0000e-04	5.00000000	True
m_SU_SG_SU_3_0_0_1	1.10915311	1.109153108617096	1.0000e-04	5.00000000	True
p_PU_SG_PU_1_n1_1_1	-0.86104140	-0.8610414024232179	-3.14159265	3.14159265	False
p_PU_SG_PU_1_1_n1_1	-0.86104140	-0.8610414024232179	-3.14159265	3.14159265	False	p_PU_SG_PU_1_n1_1_1
p_PU_SG_PU_3_n1_1_1	-3.12044446	-3.1204444620772467	-3.14159265	3.14159265	True
p_PU_SG_PU_3_1_n1_1	-3.12044446	-3.1204444620772467	-3.14159265	3.14159265	False	p_PU_SG_PU_3_n1_1_1
p_SU_SG_SU_1_0_0_1	2.61122920	2.611229196458127	-3.14159265	3.14159265	True
p_SU_SG_SU_3_0_0_1	-0.07867828	-0.07867827542158025	-3.14159265	3.14159265	True

This sets self.params from the matrix elements, which are a set of (real) parameters for lmfit, as a Parameters object.

Note that:

The input matrix elements are converted to magnitude-phase form, hence there are twice the number as the input array, and labelled m or p accordingly, along with a name based on the full set of QNs/indexes set.
One phase is set to vary=False, which defines a reference phase. This defaults to the first phase item.
Min and max values are defined, by default the ranges are 1e-4<mag<5, -pi<phase<pi.
Relationships between the parameters are set by default, but can be set manually, see section below, or pass paramsCons=None to skip.

7.4.4. Running fits#

7.4.4.1. Single fit#

With the parameters and data set, just call self.fit()!

Statistics and outputs are handled by lmfit, which includes uncertainty estimates and correlations in the fitted parameters.

# data.randomizeParams()  # Randomize input parameters if desired
                          # For method testing using known initial params is also useful
data.fit()

# Check fit outputs - self.result shows results from the last fit
data.result

Show code cell output Hide code cell output

Fit Statistics

fitting method	leastsq
# function evals	1320
# data points	156
# variables	7
chi-square	0.01615077
reduced chi-square	1.0839e-04
Akaike info crit.	-1417.40039
Bayesian info crit.	-1396.05140

Variables

name	value	standard error	relative error	initial value	min	max	vary	expression
m_PU_SG_PU_1_n1_1_1	1.81863008	0.05116438	(2.81%)	1.7846157536101068	1.0000e-04	5.00000000	True
m_PU_SG_PU_1_1_n1_1	1.81863008	0.05116438	(2.81%)	1.7846157536101068	1.0000e-04	5.00000000	False	m_PU_SG_PU_1_n1_1_1
m_PU_SG_PU_3_n1_1_1	0.78365089	0.12317226	(15.72%)	0.802904951323892	1.0000e-04	5.00000000	True
m_PU_SG_PU_3_1_n1_1	0.78365089	0.12317225	(15.72%)	0.802904951323892	1.0000e-04	5.00000000	False	m_PU_SG_PU_3_n1_1_1
m_SU_SG_SU_1_0_0_1	2.67947346	0.06072913	(2.27%)	2.686062120382649	1.0000e-04	5.00000000	True
m_SU_SG_SU_3_0_0_1	1.15266253	0.14480693	(12.56%)	1.109153108617096	1.0000e-04	5.00000000	True
p_PU_SG_PU_1_n1_1_1	-0.86104140	0.00000000	(0.00%)	-0.8610414024232179	-3.14159265	3.14159265	False
p_PU_SG_PU_1_1_n1_1	-0.86104140	0.00000000	(0.00%)	-0.8610414024232179	-3.14159265	3.14159265	False	p_PU_SG_PU_1_n1_1_1
p_PU_SG_PU_3_n1_1_1	-3.14159265	0.26152765	(8.32%)	-3.1204444620772467	-3.14159265	3.14159265	True
p_PU_SG_PU_3_1_n1_1	-3.14159265	0.26152765	(8.32%)	-3.1204444620772467	-3.14159265	3.14159265	False	p_PU_SG_PU_3_n1_1_1
p_SU_SG_SU_1_0_0_1	2.66958866	0.13105172	(4.91%)	2.611229196458127	-3.14159265	3.14159265	True
p_SU_SG_SU_3_0_0_1	0.16332776	0.09525204	(58.32%)	-0.07867827542158025	-3.14159265	3.14159265	True

Correlations (unreported correlations are < 0.100)

m_PU_SG_PU_1_n1_1_1	p_SU_SG_SU_1_0_0_1	+0.9725
m_PU_SG_PU_1_n1_1_1	m_PU_SG_PU_3_n1_1_1	-0.9647
m_SU_SG_SU_1_0_0_1	m_SU_SG_SU_3_0_0_1	-0.9630
m_PU_SG_PU_3_n1_1_1	p_SU_SG_SU_1_0_0_1	-0.9433
m_SU_SG_SU_1_0_0_1	p_SU_SG_SU_1_0_0_1	+0.8318
m_SU_SG_SU_3_0_0_1	p_SU_SG_SU_1_0_0_1	-0.8029
p_PU_SG_PU_3_n1_1_1	p_SU_SG_SU_3_0_0_1	+0.7354
m_PU_SG_PU_1_n1_1_1	m_SU_SG_SU_1_0_0_1	+0.6905
m_PU_SG_PU_3_n1_1_1	m_SU_SG_SU_1_0_0_1	-0.6741
m_PU_SG_PU_1_n1_1_1	m_SU_SG_SU_3_0_0_1	-0.6705
m_PU_SG_PU_3_n1_1_1	m_SU_SG_SU_3_0_0_1	+0.5872
m_PU_SG_PU_3_n1_1_1	p_PU_SG_PU_3_n1_1_1	-0.5076
p_PU_SG_PU_3_n1_1_1	p_SU_SG_SU_1_0_0_1	+0.4787
m_PU_SG_PU_1_n1_1_1	p_PU_SG_PU_3_n1_1_1	+0.4280
m_PU_SG_PU_3_n1_1_1	p_SU_SG_SU_3_0_0_1	-0.4227
m_SU_SG_SU_1_0_0_1	p_PU_SG_PU_3_n1_1_1	+0.3819
m_PU_SG_PU_1_n1_1_1	p_SU_SG_SU_3_0_0_1	+0.3771
p_SU_SG_SU_1_0_0_1	p_SU_SG_SU_3_0_0_1	+0.3208
m_SU_SG_SU_3_0_0_1	p_PU_SG_PU_3_n1_1_1	-0.2783

# Plot results with data
# data.BLMfitPlot(backend='hv')
data.BLMfitPlot()

Dataset: subset, AFBLM
Dataset: 0, AFBLM

../_images/a504585a66aa694e2981684ead7c70a31555fec75ec5d05a6ffe3a2266f35693.png

7.4.4.2. Extended execution methods, including parallel and batched execution#

See https://pemtk.readthedocs.io/en/latest/fitting/PEMtk_fitting_demo_multi-fit_tests_130621-para_010922.html

(1) serial execution

Either:

Manually with a loop.
With self.multiFit() method, although this is optimised for parallel execution (see below).

import time

start = time.time()

# Maual execution
for n in range(0,10):
    data.randomizeParams()
    data.fit()
    
end = time.time()
print((end - start)/60)
    
# Or run with self.multiFit(parallel = False)
# data.multiFit(nRange = [0,100], parallel = False)

0.7577979763348898

# We now have 10 fit results
data.data.keys()

dict_keys(['subset', 'orb6', 'orb5', 'ADM', 'pol', 'sim', 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10])

(b) parallel execution

Updated version including parallel fitting routine with self.multiFit() method.

This currently uses the XYZpy library for quick parallelization, although there is some additional setup overhead in the currently implementation due to class init per fit batch. The default aims to set ~90% CPU usage, based on core-count.

# Multifit wrapper with range of fits specified
# Set 'num_workers' to override the default.
data.multiFit(nRange = [0,10], num_workers=20)

Number of processors:  64 
Running pool on:  20
* sparse not found, sparse matrix forms not available. 
* natsort not found, some sorting functions not available. 
* Setting plotter defaults with epsproc.basicPlotters.setPlotters(). Run directly to modify, or change options in local env.
* Set Holoviews with bokeh.
* pyevtk not found, VTK export not available. 
* sparse not found, sparse matrix forms not available. 
* natsort not found, some sorting functions not available. 
* Setting plotter defaults with epsproc.basicPlotters.setPlotters(). Run directly to modify, or change options in local env.
* Set Holoviews with bokeh.
* pyevtk not found, VTK export not available. 
* sparse not found, sparse matrix forms not available. 
* natsort not found, some sorting functions not available. 
* Setting plotter defaults with epsproc.basicPlotters.setPlotters(). Run directly to modify, or change options in local env.
* Set Holoviews with bokeh.
* pyevtk not found, VTK export not available. 
* sparse not found, sparse matrix forms not available. 
* natsort not found, some sorting functions not available. 
* Setting plotter defaults with epsproc.basicPlotters.setPlotters(). Run directly to modify, or change options in local env.
* Set Holoviews with bokeh.
* pyevtk not found, VTK export not available. 
* sparse not found, sparse matrix forms not available. 
* natsort not found, some sorting functions not available. 
* Setting plotter defaults with epsproc.basicPlotters.setPlotters(). Run directly to modify, or change options in local env.
* Set Holoviews with bokeh.
* pyevtk not found, VTK export not available. 
* sparse not found, sparse matrix forms not available. 
* natsort not found, some sorting functions not available. 
* Setting plotter defaults with epsproc.basicPlotters.setPlotters(). Run directly to modify, or change options in local env.
* Set Holoviews with bokeh.
* pyevtk not found, VTK export not available. 
* sparse not found, sparse matrix forms not available. 
* natsort not found, some sorting functions not available. 
* Setting plotter defaults with epsproc.basicPlotters.setPlotters(). Run directly to modify, or change options in local env.
* Set Holoviews with bokeh.
* pyevtk not found, VTK export not available. 
* sparse not found, sparse matrix forms not available. 
* natsort not found, some sorting functions not available. 
* Setting plotter defaults with epsproc.basicPlotters.setPlotters(). Run directly to modify, or change options in local env.
* Set Holoviews with bokeh.
* pyevtk not found, VTK export not available. 

  0%|                                                                                                              | 0/10 [00:00<?, ?it/s]OMP: Info #271: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead.
OMP: Info #271: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead.
OMP: Info #271: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead.
OMP: Info #271: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead.
OMP: Info #271: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead.
OMP: Info #271: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead.
OMP: Info #271: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead.
OMP: Info #271: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead.
OMP: Info #271: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead.
OMP: Info #271: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead.
100%|#####################################################################################################| 10/10 [00:46<00:00,  4.62s/it]

(c) Dump data

Various options are available. The most complete is to use Pickle (default case), although this is not suggested for archival use. For details see https://epsproc.readthedocs.io/en/dev/dataStructures/ePSproc_dataStructures_demo_070622.html

outStem = 'dataDump_N2'  # Set for file save later

# data.writeFitData(fName='N2_datadump')  # Use 'fName' to supply a filename
data.writeFitData(dataPath = dataPath, outStem=outStem)  # Use 'outStem' to define a filename which will be appended with a timestamp
                                                        # Set dataPath if desired, otherwise will use working dir

Dumped self.data to /home/jovyan/jake-home/buildTmp/_latest_build/html/doc-source/part2/n2fitting/dataDump_N2_280423_09-40-05.pickle with pickle.
Dumped data to /home/jovyan/jake-home/buildTmp/_latest_build/html/doc-source/part2/n2fitting/dataDump_N2_280423_09-40-05.pickle with pickle.

PosixPath('/home/jovyan/jake-home/buildTmp/_latest_build/html/doc-source/part2/n2fitting/dataDump_N2_280423_09-40-05.pickle')

7.5. Quick setup with script#

The steps demonstrated above are also wrapped in a helper script, although some steps may need to be re-run to change selection properties or ranges.

# Run general config script with dataPath set above
%run {dataPath/"setup_fit_demo.py"} -d {dataPath}